A discrete statistical mechanics approach to aeolian ripple dynamics 
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The interaction between a fluid and a granular material plays a crucial role in a large class of 
phenomena such as landscape morphology and transport of sediments, aeolian sand dunes formation 
and ripples dynamics. Standard models involve deterministic continuum equations or, alternatively, 
Lattice Boltzmann and Lattice Gas Cellular Automata. We here introduce a toy-model to address 
the issue of aeolian ripple formation and evolution. Our simplified approach accounts for the basic 
physical mechanisms and enables to reproduce the observed phenomenology in the framework of an 
innovative statistical mechanics formulation. 
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INTRODUCTION AND GENERAL 
BACKGROUND 



The complex interactions between a granular mate- 
rial and a fluid may lead to the formation of spectacular 
structures, like erosion patterns, dunes and ripples. We 
are here interested in ripple patterns that emerge due to 
the interactions between a fluid (air) flow and a granu- 
lar material (sand). The conditions of the flow and the 
characteristics of the granular material give origin to a 
large variety of morphologies 0, Q . 

Ripples have been recently photographed even on 
Mars, qualitatively resembling the ones observed on the 
Earth surface, see Fig.Q] From a direct inspection of the 
peculiar shapes of such coherent structures, one could in 
principle aim at reconstructing the characteristics of the 
surrounding environment. It is therefore of paramount 
importance to develop dedicated models enabling to suc- 
cessfully address the issue of ripple formation and evo- 
lution. Phenomcnologically, the correct interpretative 
framework is provided by the pioneering investigations 
by Bagnold [l| based on wide observational and experi- 
mental facts. In particular, the two key processes, namely 
saltation and surface creeping (reptation), are outlined in 
Rcf. 0. 

Saltation refers to the process by which a grain is en- 
trained by the fluid, accelerated and transported to an- 
other site. In wind flows the transition from rest to en- 
trainment occurs as a sudden jump, when the lift force 
is able to dislodge a particle from its gravitational trap. 
Due to the large density differences between sand and air, 
this occurs for large shear velocities. The takeoff angle 
(with respect to horizontal) is large, and so is the flight 
time and transport. The subsequent impact with the 
other particles occurs at relatively high energy, and even 
when the restitution coefficient is small, this energy may 
be sufficient to dislodge other particles. These reptated 
particles are then transported by the fluid, but their ini- 
tial velocity is generally lower than that associated to the 
saltating particles, thus preventing avalanche effects. 
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FIG. 1: (color online) Pattern of aeolian ripples in south 
central Colorado (upper) 0. Giant sand dunes and ripples, 
steeper than on Earth, has been recently observed on Mars 
(NASA) (lower) 0]. 



Reptated particles finally stop in a local gravitational 
trap, loosing energy via anelastic collisions. One of the 
standard assumption [j| is that this reptation mechanism 
is the main origin of the creeping flow, including what in 
other context would be referred as toppling. 

Which one of these phenomena is more relevant, de- 
pends on the combined effects of drag and lift forces, 
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gravity and dissipation of energy during collisions. How- 
ever the main difficulty arises in modeling the interaction 
between the fluid and the grains. This interaction de- 
pends on the distribution of the fluid velocity, that may 
vary widely in space and time, on the shape of the grains 
and on their relative density, that gives origin to buoy- 
ancy. Moreover, the instantaneous fluid velocity field de- 
pends on the profile of the sand bed, which is modified 
by erosion and deposition, and on the shielding effect of 
saltating particles. Self-organizedpatterns are originated 
by these nontrivial interactions [y, Q ■ Dedicated models 
aim at reproducing some of the most relevant aspects 
of these phenomena. In particular, following Ref. 0, 
ripples are characterized by being asymmetric and their 
formation and persistence is determined by flux intensity. 
Small ripples are observed to travel faster than large ones 
and the temporal evolution of the maximum ripple height 
is limited and not linear. 

In order to study the erosion/deposition process con- 
tinuum models of mechanics have been proposed 0, 0, 
ITnj . In addition accurate hydrodynamical descriptions 
have been studied by various authors 0, O, 0, El- 
These alternative approaches enable to quantitavely 
study the process of ripple formation by incorporating a 
detailed representation of the fluid flow via Navier-Stokcs 
equation, and/or facing the problem of considering the 
flow near the soil. However, the latter models arc com- 
putationally expensive and difficult to treat analitically 
particularly in presence of complex and time-dependent 
boundaries, or when turbulence and other fluctuating as- 
pects play an active role. However, linear stability anal- 
ysis [(j B EJ enables to predict the instability regime of 
a flat surface and quantifies its associated growth rate. 

To gain more insight into the crucial interplay between 
erosion and deposition, beyond the linear approximation, 
a series of simplified theoretical frameworks (toy models) 
have been developed, focusing only on those aspects sup- 
posed to be the relevant ones. Following these lines, it is 
customary to replace the complex fluid velocity distribu- 
tion with a limited number of aggregated data, like the 
average shear velocity and the associated fluctuations. 

Simple models of sand ripples dynamics were first 
introduced by Anderson 0, UJ for grain segregation 
and stratigraphy. A discrete stochastic model was fur- 
ther proposed about a decade ago by Werner and Gille- 
spie [lfj. Another minimal model, widely adopted in the 
relevant literature, was proposed the same year by Nishi- 
mori and Ouchi The latter, termed NO, qualitatively 
captures some of the peculiar features of ripple evolu- 
tion. Within this scenario, the saltation and reptation 
are accounted for and shown to produce the spontaneous 
formation of a characteristic ripple patterns. Though it 
represents a significant step forward in the comprehen- 
sion of the basic mechanism underlying the phenomenon, 
the NO model allows for non realistic structures of in- 
finite heights, since in this model there is no explicit 
or implicit mechanism that leads to the appearance of 
a critical angle of repose. Recently, Caps and Vandc- 



wallc 17] modified the preexisting scheme by including 
explicitly the effect of avalanches (SCA model: Salta- 
tion Creep and Avalanches). This modification results in 
asymmetric ripple profiles and induces a saturation for 
the maximal height. 

Cellular automata modeling of sand transportation 
was introduced by Anderson and Bunnan |l8j and by 
Werner and Gillespie 0|. In these models, the driving 
mechanism for sand transport is the saltation/reptation 
dynamics, eventually complemented by toppling, that 
corresponds to diffusion in a continuous model. Masselot 
and Chopard 0| also introduced a cellular automata for 
snow and sand transportation. They explicitly modeled 
the fluid flow by means of lattice Boltzmann methods, 
while the granular phase is represented as a probabilis- 
tic cellular automaton. The erosion mechanism here is 
modeled by a constant probability of detachment, and 
local rearrangements are again achieved by a toppling 
mechanism. 

The elementary building blocks of these stochastic 
models, like the erosion, deposition and toppling steps, 
have a phenomcnological nature, implying that the prob- 
ability of their occurrence has to be measured experi- 
mentally. In this work we propose to adopt a standard 
statistical mechanics description of the elementary steps, 
coupled with an external forcing that drives the system 
out of equilibrium. The role of temperature is here played 
by the fluctuations of the velocity field. 

We focus on the formation and evolution of aeolian rip- 
ples. The proposed local-equilibrium dynamics enables 
to reproduce the main characteristics of ripple dynam- 
ics, like the observed stable states and the saturation of 
ripples heights, without including an explicit mechanism 
for toppling or other local rearrangements. This simple 
scheme may be easily extended to include a more detailed 
description of a real fluid. 

The paper is organized as follows: in section [H] the 
model is introduced. Section IIIII is devoted to discuss 
the numerical implementation. Results are presented in 
Section llVl Finally in Section[V]we draw our conclusion. 



II. THE MODEL 

We consider a one dimensional discrete model, L be- 
ing the extension of the segment partitioned in N equally 
spaced intervals, and assuming periodic boundary condi- 
tions, Fig.H 

To each site i we associate the height hi which labels 
the number of particles constituting the i th slice of the 
sand bed. Further, we consider a bunch of rii particles 
flowing over the bed. The system is therefore composed 
of two interacting layers of particles. 

The two processes driving the dynamics of the system 
are: erosion, which occurs when a resting grain belong- 
ing to the surface of the sediments layer is entrained by 
the fluid and deposition, that mimics the deposition of a 
flowing grain. 
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FIG. 2: Sand bed schematization in the ID model. For a 
given site i of height hi, the symmetric interval of neighbors 
of amplitude R is shown (dotted line). In order to account 
for flux asymmetry the interval considered is shifted of an 
amount S (dashed line). 
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FIG. 3: In order to account for bed geometry and hydro- 
dynamics effects, we suppose a grain on the bed surface to 
experience a force proportional to the difference between its 
height and the mean height of a local interval of neighbors. 
In particular, the erosion (deposition) probability increases as 
the height of the site being considered is higher (lower) than 
the mean height of the neighboring sites. 



energy scales therefore as: 

Ei = (h. 



(1) 



The erosion and deposition processes are hence charac- 
terized by means of the following change in energy: 



Erosion 


Deposition 


h t (t + l) = hi(t) - 1 
AE erA = 2{h m -fn + 0.5) 


hi(t + l) = hi(t) + 1 
AE dep>i = 2{fu -h m + 0.5) 



Consequently, it is reasonable to assume the erosion 
and deposition probabilities pi) : 



P 



dep.? 



1 if AE ei>i < 0, 

exp(— [3 e AE er ^) otherwise; 

'l if A£ dep ,i < 0, 

exp(— /S^Ai^dop.i) otherwise. 



where (3 e and (3d are constant parameters, analogous to 
the inverse of effective temperatures 1/T e , 1/T d . The 
system is largely out of equilibrium, and these tempera- 
tures are assumed to be related to the amplitude of the 
fluctuations of the velocity field in correspondence of the 
typical erosion and deposition events. 



III. NUMERICAL IMPLEMENTATION 

The evolution rule, applied simultaneously to each site, 
is based on a Metropolis Monte Carlo dynamics |2l| : first 
the deposition is made to happen followed by the subse- 
quent erosion step. Focusing on the i th site, the deposi- 
tion step yields: 



The evolution of the system is then modeled as fol- 
lows: focusing on the i th site, we select an interval of R 
neighbors, asymmetrically shifted by an amount S, see 
Fig. |21 The amplitude of the interval accounts for the 
range of local interactions and is shown to be correlated 
to the shape of the ripple. The quantity S is introduced 
to model hydrodynamics effects, such as lift and drag 
forces 0, 123, which result in an asymmetry of the 
flow. 

We then calculate the mean height of the selected in- 
terval, h m . As a reasonable hypothesis, assume that the 
probability to experience an erosion event increases by 
augmenting the gap between hi and h m , under the con- 
straint hi > h rn . Conversely, the deposition will most 
probably occur when the positive difference h m — hi gets 
larger, Fig. [21 

Following the previous reasoning, as a first approxima- 
tion, we suppose the force acting on a particle to depend 
on the difference between the height of the site being 
considered and the mean height of the neighbors. The 



1 . for each of the rii flowing particles a uniformly dis- 
tributed random number r is extracted; 

2. if r < Pd deposition occurs and the height of the 
site is increased by one, hi(t + 1) = hi(t) + 1, while 
the bunch of flowing particles is decreased by one, 
rii(t+ 1) = rii{t) - 1, 

The erosion is characterized by: 

1 . a uniformly distributed in the unit interval random 
number r is generated; 

2. according to the Monte Carlo Metropolis rule the 
site is eroded if r < P e ; in this case the height of 
the site is decreased by one, hi(t + 1) = hi(t) — 1, 
and the pool of eroded particle is increased by one, 
rii{t+ 1) = m(t) + 1. 

The procedure is iterated and the evolution of the 
heights monitored. 
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FIG. 5: When two ripples of similar sizes encounter, an ex- 
change in their relative positions occur (take into account 
periodic boundary conditions). 
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mapshots of ripple dynamics. Note that the 
axes are different. 



IV. RESULTS 

Numerical simulations are performed starting from an 
initial uniformly random generated river-bed. Small in- 
homogeneities are enhanced as time progresses, and even- 
tually result in macroscopic ripples that display a char- 
acteristic asymmetric profile. 

A sequence of successive snapshots of the dynamics is 
collected in Figure 01 and allows to qualitatively investi- 
gate the process of formation of coherent structures. 

The displacement of a ripple is a consequence of the 
combined effects of erosion and deposition: the grains 
are eroded in the stoss side and deposed in the lee one. 
This is a dynamical mechanism: there is a continuous 
exchange between the rest particles and the flowing ones. 
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FIG. 6: An encounter between ripples of different sizes results 
in the merging of the structures. 



The global consequence of this flow is that the lower end 
of the stoss side is eroded and this matter deposes in the 
lee part, leading to the displacement of the ripple with a 
velocity that decreases with size. 

The interaction mechanism of two ripples is rather 
complex, as illustrated by the "collision" of two ripples 
of similar size, Fig- EI I n order to collide, the front ripple 
has to be larger than the rear one, which is consequently 
faster. When the two ripples approach, the lower part of 
the stoss side of the front ripple is no more eroded. This 
is due (in the model) to the increasing average height, 
that includes the contribution of the approaching rear 
ripple, fn the reality, this corresponds to the reduction 
in erosion due to the shielding effect of the rear ripple. 

The region of the stoss side next to this one becomes 
the source of eroded particles, thus forming a local sink. 
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FIG. 7: Ripple schematization, where xi (x s ) is the projection 
of the lee (stoss) slope and h is the ripple height. 



The size of the front ripple is decreased and its speed 
increased. If the relative difference of velocity of the two 
ripples is low enough (i.e. for similar ripples of similar 
size) this depression may proceed enough to make the 
front ripple detach from the rear one, Fig. [3J 

If the rear ripple is small, the sinking region is con- 
tinuously moved downwind and finally the two ripples 
coalesce. This mechanism is illustrated in Fig. [S] In the 
third panel (t = 825) one can still recognize a signature 
of this process: focusing on the right ripple, originated 
by the interaction of two ripples of different size (see first 
panel) , one can still identify the protruding bump on the 
lee side. This is the relic of the highest peak displayed 
by the right ripple in the second panel, that experienced 
a reduction in size and consequently proceeded faster. 
However, this bump is eventually screened by the rear 
portion of the ripple. As a consequence it stops and is 
therefore engulfed in the incoming massive bulk. 

These observations are in agreement with direct mea- 
surements and provide a first validation of our simplified 
interpretative framework []]. 

Label with I the linear size of a typical ripple and as- 
sume h to represent its characteristic height. By tun- 
ing the parameter R, i.e. varying the extension of the 
segment that defines the interacting region, one mod- 
ulates the ratio h/l and operates an a priori selection 
among various types of structures (ripples, megaripplcs, 
giantripples) based on their intrinsic geometry. The cru- 
cial role of R has been investigated through a dedicated 
campaign of simulations: by assigning larger values of 
R corresponds to generate less peaked structure, which 
translates into a systematic reduction of the quantitative 
indicator h/l. Focusing on ripples, we assume R = 32, 
and speculate on the role of the remaining parameters S 
and (3d in the formation of the ripple |22l | : few quantities 
of paramount importance are monitored and compared 
with analogous predictions reported in the classical liter- 
ature HE3. 

As already anticipated, within our simplified scheme, 
the ripples present an asymmetric shape which can be 
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FIG. 8: Time evolution of ripple asymmetry a as a function 
of the ratio S/R. Left-inset: asymptotic value of ripple asym- 
metry as a function of S/R. Right-inset: plot of the dynamic 
angle of repose 8 r against the ratio S/R (circles: numerical 
values; solid line: numerical fit). 
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FIG. 9: Density of particles p as a function of the deposition 
distance, for (3d — 0.1 and 1.0. Inset: characteristic length 
scale A of deposition as a function of (3d- 



measured by introducing the aspect ratio a: 



where x s (resp. xi) stands for the projection of the stoss 
(resp. lee) slope on the horizontal axis, as depicted in 
Figure [7| If a = 1 the ripples are symmetric, while a / 
implies an asymmetry. In the main panel of Figure |SJ 
the evolution of the ripple aspect ratio a is plotted as 
function of time for different values of the ratio S/R. 
An initial growth is observed, followed by a subsequent 
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FIG. 10: Exponentially increase of the maximum ripple height 
hmax with time for different values of the ratio S/R. 



saturation towards an asymptotic plateau, o~ as . A self- 
consistent selection mechanism is therefore operated by 
the system and eventually only one specific class of rip- 
ples arises and occupies the one dimensional lattice. The 
following ansatz, 

a = a as (1 - cxp(-at)) , (3) 

is numerically fitted to the simulated profiles of Figure |S] 
and shown to interpret well the data. To better visual- 
ize the tendency of enhancing the degree of asymmetry 
for increased values of the local distortion S (working at 
constant R), a as is represented versus the ratio S/R in 
the top- left inset of Figure |SJ Further, to provide a com- 
plete characterization of the morphology of the ripple, 
we calculate the repose angle 8 r , dynamically selected 
within our proposed approach, as function of the control 
quantity S/R. 

The temporal evolution of a was previously monitored 
in Ref. [13 for both the NO and SCA models. The orig- 
inal NO formulation predicts almost symmetric profiles 
and therefore a ~ 1. Conversely, for the case of the SCA 
o~ grows linearly in time and then relaxes to a final value. 
This remarkable improvement was achieved by Caps and 
Vandewalle by postulating the existence of a repose angle 
r and modeling the process of avalanches, not included 
in the NO philosophy. It is worth emphasizing that a 
somehow similar mechanism (the saturation for a being 
here exponentially approached) is here reproduced with- 
out invoking a priori the existence of a limiting angle. 

In Figure the deposition process is investigated to 
shed light into the crucial role of (3d- The normalized 
density of particles p as function of the deposition dis- 
tance is plotted for distinct values of (3d- The measured 
distributions decays exponentially and the characteristic 
length scale A is represented versus (3d in the small inset. 
As expected, for larger values of (3d the particles spend 
more time in the surrounding halo and retard the depo- 
sition event. We therefore suggest that (3d provides an 
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indirect control of the characteristic length of the repta- 
tion process. 

Finally, we studied the dynamical evolution of the 
maximum ripple height h max . In |l7j the SCA model was 
shown to reproduce the non linear evolution of the ripple 
amplitude h ma x, this success being ascribed to the new 
ingredients introduced with respect to the NO scenario. 
Results of our simulations arc reported in Figure ll(JI as 
for the SCA an exponential growth law is also found in 
the framework of our probabilistic approach, thus rein- 
forcing its validity as an alternative tool to address the 
relevant issue of ripple formation. 



produce important observed features, despite its intrinsic 
simplicity. In particular, the irreversible and not triv- 
ial coarsening dynamics with merging and scattering of 
structures, the saturating value of the maximum height 
of the ripples and the asymmetry of ripple structures 
have been reproduced. The results are critically com- 
pared with the classical literature 0, 0] , outlining the 
role of the parameters involved in our formulation and 
their physical interpretation. This simple formulation 
may be easily extended to include a more detailed de- 
scription of the fluid flow and the characteristics of the 
granular phase. 



V. CONCLUSION 



The proposed local-equilibrium model for the study 
of aeolian ripple dynamics has shown to successfully re- 
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